#define VERBOSE_GAUSS_ELIMINATION

using System;
using System.Collections.Generic;
using System.Text;


namespace Algorithms.Numerical
{


    public class Algebra
    {
        public static Algebra Instance
        {
            get
            {
                if (instance == null)
                    instance = new Algebra();
                return instance;
            }
        }

        /**
         *  Calculate base matrix of solutions of:
         *      MX=0
         *  Columns of resultant matrix form a base of
         *  solutions space. This matrix can be degenerated
         *  and have width of 0. In such case solution space
         *  is vector 0.
         * 
         * 
         */
        public Matrix SolveUniform( Matrix coefs, bool inPlace )
        {
            int i;  // general purpose counter

            GaussMatrix m = GaussEliminateBelow(coefs, inPlace);
            m = GaussEliminateAbove(m, true);
            #if VERBOSE_GAUSS_ELIMINATION
            Console.WriteLine("After elim. above:");
            Console.Write(m.ToString());
            #endif

            #region Create selection arrays
            int[] indepColsIdx = new int[m.Rank];
            int[] depColsIdx = new int[m.Width - m.Rank];
            int last = -1;
            int j = 0;
            for (i = 0; i < m.Rank; i++)
            {
                for (int k = last + 1; k < m.IndepColIdx(i); k++)
                    depColsIdx[j++] = k;
                last = indepColsIdx[i] = m.IndepColIdx(i);
            }
            for (int k = last + 1; k < m.Width; k++)
                depColsIdx[j++] = k;
            #endregion

            Matrix indep = m.SelectionView(null, indepColsIdx);
            Matrix dep = m.SelectionView(null, depColsIdx).Copy();

            #if VERBOSE_GAUSS_ELIMINATION
            Console.WriteLine("Indep:");
            Console.Write(indep.ToString());

            Console.WriteLine("Dep:");
            Console.Write(dep.ToString());
            #endif

            m = MakeIndepColumnsUniform(m);    

            #if VERBOSE_GAUSS_ELIMINATION
            Console.WriteLine("Indep after uniformization:");
            Console.Write(indep.ToString());
            #endif
        
            
            Matrix r = MatrixFactory.Instance.MakeZero(m.Width, m.Width - m.Rank);
            r.SelectionView(indepColsIdx, null).Assign(m.SubmatrixView(0, 0, m.Rank, m.Width)
                .SelectionView(null, depColsIdx));
            r.Mul(-1);
            r.SelectionView(depColsIdx, null).Assign(MatrixFactory.Instance.Mult(indep[0,0],
                MatrixFactory.Instance.MakeIdentity(m.Width - m.Rank, m.Width - m.Rank), true));

        
            #if VERBOSE_GAUSS_ELIMINATION
            Console.WriteLine("Res before:");
            Console.Write(r.ToString());
            #endif

            ReduceValuesInRows(r.TransposedView);

            return r;
        }


        public GaussMatrix GaussEliminateBelow(Matrix m, bool inPlace)
        {
            if ( !inPlace)
                m = m.Copy();

            int[] auxIndepIdx = new int[m.Width];
            int auxIndepCnt = 0;

            // For each column
            int x = 0;
            int y = 0;
            while (x < m.Width && y < m.Height)
            {
                // In column x find first non-zero element at or below position y
                int i = y;
                while (i < m.Height && m[i, x] == 0)
                    i++;

                // If found
                if (i < m.Height)
                {
                    // This will be one of lin. indep. columns. Remember it
                    auxIndepIdx[auxIndepCnt++] = x;

                    // If element found not at diagonal
                    if (i > y)
                    {
                        // move it to diagonal
                        m = m.ChangedRowsView(i, y);
                        #if VERBOSE_GAUSS_ELIMINATION
                        Console.WriteLine("Line changed");
                        Console.Write(m.ToString());
                        #endif
                    }

                    #if VERBOSE_GAUSS_ELIMINATION
                    Console.WriteLine(string.Format("X: {0}; Y: {1}, I: {2}", x, y, i));
                    #endif

                    // Eliminate elements below diagonal
                    Matrix y_row = m.SelectionView(new int[] { y }, null);
                    for (i = y + 1; i < m.Height; i++)
                    {
                        double y_mult = -m[i, x] / m[y, x];
                        double i_mult = 1;
                        // Try to keep integers if posible
                        if (IsInt(m[i, x]) && IsInt(m[y, x]))
                        {
                            int gcd = Gcd((int)m[i, x], (int)m[y, x]);
                            y_mult = -m[i, x] / gcd;
                            i_mult = m[y, x] / gcd;
                        }

                        // M_i = i_mult * M_i + y_mult * M_y
                        Matrix i_row = m.SelectionView(new int[] { i }, null);
                        i_row.Mul(i_mult);
                        i_row.Add(MatrixFactory.Instance.Mult(y_mult, y_row, false));
                    }

                    //!!
                    y++;
                } // else
                // Process next column
                x++;
                #if VERBOSE_GAUSS_ELIMINATION
                Console.WriteLine("Step done");
                Console.Write(m.ToString());
                #endif
            }

            int[] indepIdx = new int[auxIndepCnt];
            for (int k = 0; k < auxIndepCnt; k++)
                indepIdx[k] = auxIndepIdx[k];

            return new GaussMatrix(ReduceValuesInRows(m), indepIdx);
        }

        public GaussMatrix GaussEliminateAbove(GaussMatrix m, bool inPlace)
        {
            if (!inPlace) {
                int[] aux = new int[m.Rank];
                for(int auxi=0; auxi<m.Rank; auxi++ )
                    aux[auxi] = m.IndepColIdx(auxi);
                m = new GaussMatrix(m.Copy(), aux);
            }

            for (int y = 0; y < m.Rank; y++)
            {
                int x = m.IndepColIdx(y);

                // eliminate all above
                Matrix y_row = m.SelectionView(new int[] { y }, null);
                for (int i = y - 1; i >= 0; i--)
                {
                    double y_mult = -m[i, x] / m[y, x];
                    double i_mult = 1;
                    // Try to keep integers if posible
                    if (IsInt(m[i, x]) && IsInt(m[y, x]))
                    {
                        int gcd = Gcd((int)m[i, x], (int)m[y, x]);
                        y_mult = -m[i, x] / gcd * Math.Sign(m[y, x]);
                        i_mult = Math.Abs(m[y, x]) / gcd;
                    }
                    // M_i = i_mult * M_i + y_mult * M_y
                    Matrix i_row = m.SelectionView(new int[] { i }, null);
                    i_row.Mul(i_mult);
                    i_row.Add(MatrixFactory.Instance.Mult(y_mult, y_row, false));
                }
            }
            ReduceValuesInRows(m);
            return m;
        }


        /**
         * If values in each row are integers and has common divisior,
         * they are reduced.
         * Works in place.
         */
        public Matrix ReduceValuesInRows(Matrix m)
        {
            for (int y = 0; y < m.Height; y++)
            {
                int gcd = -1;
                int x;
                for ( x = 0; x < m.Width; x++ )
                    if (IsInt(m[y, x]))
                    {
                        if (gcd == -1)
                            gcd = (int)Math.Abs(m[y, x]);
                        else
                            gcd = Gcd(gcd, (int)Math.Abs(m[y, x]));
                    }
                    else
                        break;
                // All values in row are integers
                if (x == m.Width && gcd > 0)
                    m.SelectionView(new int[] { y }, null).Mul(1.0 / gcd);
            }
            return m;
        }

        public GaussMatrix MakeIndepColumnsUniform(GaussMatrix m)
        {
            int i;
            int gcd = -1;
            int lcm = 1;
            for (i = 0; i < m.Rank; i++)
            {
                if ( IsInt( m[i, m.IndepColIdx(i)]))
                {
                    if (i == 0)
                        gcd = (int)m[i, m.IndepColIdx(i)];
                    else
                        gcd = Gcd( gcd, (int)Math.Abs( m[i, m.IndepColIdx(i)] ));
                    lcm *= (int)Math.Abs(m[i, m.IndepColIdx(i)]);
                } else
                    break;
            }

            // All integers
            double numerator = 1;
            if (i == m.Rank)
            {
                lcm /= gcd;
                numerator = lcm;
            }

            for (i = 0; i < m.Rank; i++)
            {
                double mul = numerator / m[i, m.IndepColIdx(i)];
                m.SelectionView(new int[] { i }, null).Mul(mul);
            }

            return m;
        }

        private bool IsInt(double d)
        {
            return Math.Truncate(d) == d;
        }

        private int Gcd(int a, int b)
        {
            a = Math.Abs(a);
            b = Math.Abs(b);
            while ( b != 0 )
            {
                int t = b;
                b = a % b;
                a = t;
            }
            return a;
        }


        protected Algebra()
        {
        }
        private static Algebra instance = null;
    }


    /**
     *  Matrix with additional information related to gauss elimination
     */
    public class GaussMatrix : Matrix
    {
        private Matrix m;
        private int[] indepColsIdx;

        public GaussMatrix(Matrix m, int[] indepColsIdx)
        {
            this.m = m;
            this.indepColsIdx = indepColsIdx;
        }

        public override int Height
        {
            get
            {
                return this.m.Height;
            }
        }

        public override int Width
        {
            get
            {
                return this.m.Width;
            }
        }

        public override double this[int y, int x]
        {
            get
            {
                return this.m[y, x];
            }
            set
            {
                this.m[y, x] = value;
            }
        }

        public int Rank
        {
            get
            {
                return this.indepColsIdx.Length;
            }
        }


        /*
         * Return index of n-th column from 
         * max. lin. indep. columns set.
         */
        public int IndepColIdx(int n)
        {
            return this.indepColsIdx[n];
        }

        public override string ToString()
        {
            string s = "[ \n";
            s += "  [";
            int z = 0;
            for (int x = 0; x < Width; x++)
            {
                if (z < indepColsIdx.Length && x == indepColsIdx[z])
                {
                    z++;
                    s += "****";
                }
                else
                    s += "____";
                
                
                if (x < (Width - 1))
                    s += ", ";
            }
            s += " ]\n";
            for (int y = 0; y < Height; y++)
            {
                s += "  [";
                for (int x = 0; x < Width; x++)
                {
                    s += this[y, x].ToString("N");
                    if (x < (Width - 1))
                        s += ", ";
                }
                s += " ]\n";
            }

            s += "]\n";

            return s;
        }

    }
}
